Fluctuating initial conditions in heavy-ion collisions from the Glauber approach* 
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In the framework of the Glauber approach applied to the initial stage of ultra-relativistic heavy-ion 
collisions we analyze the shape parameters of the early-formed system (fireball) and their event-by- 
event fluctuations. We test a variety of models: the conventional wounded nucleon model, a model 
admixing binary collisions to the wounded nucleons, a model with hot spots, as well as the hot-spot 
model where the deposition of energy occurs with a superimposed probability distribution. We look 
in detail at the so-called participant harmonic moments, e*, obtained by an averaging procedure 
where in each event the system is translated to its center of mass and aligned with the major 
principal axis of the ellipse of inertia. Quantitative comparisons indicate substantial relative effects 
for e* in variants of Glauber models. On the other hand, the dependence of the scaled standard 
deviation Ae*/e* on the chosen model is weak. For all models the values range from about 0.5 for the 
central collisions to about 0.3-0.4 for peripheral collisions, both for the gold-gold and copper-copper 
collisions. They are dominated by statistics and change only by 10-15% from model to model. We 
provide an approximate analytic expansion for the harmonic moments and their fluctuations given in 
terms of the fixed-axes moments. For central collisions and in the absence of correlations it gives the 
simple formula Ae* je* ~ y/A/it — 1 = 0.52. Similarly, we obtain expansions for the radial profiles of 
the higher harmonics. We investigate the relevance of the shape-fluctuation effects for jet quenching 
and find them important only for very central events. Finally, we make some comments of relevance 
for hydrodynamics, the elliptic flow and its fluctuations. We argue how smooth hydrodynamics 
leads to the known result V4 ~ v\, and further to the prediction Av4,/v4 = 2A«2/«2. 
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I. INTRODUCTION 

It was realized a few years ago in event-by-event hy- 
drodynamic studies 0, 0] of relativistic heavy-ion colli- 
sions that fluctuations of the initial shape of the fire- 
ball formed in the early stage of the reaction lead to 
quantitatively important effects for azimuthal asymme- 
try [!, 0, H, [g] . These effects, resulting from the shift of 
the center-of-mass and the rotation of the the quadrupole 
principal axis, can be seen in the analyses of the elliptic 
flow @, H, H, [HI. The purpose of this paper is to in- 
vestigate this phenomenon in detail in the framework of 
various Glauber-like approaches describing the deposi- 
tion of energy in the system in the early stages of the 
collision. Our study focuses on both the understanding 
of the statistical nature of the results, as well as on com- 
parisons of various models. The main outcome presented 
in this paper is twofold: first, we provide the Fourier 
moments and radial profiles of the so-called participant 
type, i.e. obtained with an averaging procedure where in 
each event the system is translated to its center of mass 
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and aligned with the major principal axis. Second, under 
reasonable approximations we derive analytic expansions 
which explain the basic features of the Fourier moments 
and profiles. 

The fact that the initial shape of the fireball fluctuates 
from event to event is certainly not surprising. Clearly, 
a finite number of sources consisting of wounded nucle- 
ons, binary collisions, etc., which deposit the transverse 
energy in the system at mid-rapidity do not fill the avail- 
able coordinate space uniformly due to statistical fluc- 
tuations. For instance, the center of mass of a system 
of uncorrelated particles fluctuates with a standard de- 
viation proportional to l/^/n. Similarly, the orientation 
of principal axis of the quadrupole and higher harmonic 
moments fluctuates from event to event. The statisti- 
cal component of the harmonic moments assumes aver- 
age values proportional l/^/n, where n is the number of 
sources. Since the number of sources is not so large, rang- 
ing from a few to a few hundred, these fluctuations may 
easily reach a value of a few percent or higher, large for 
studies of azimuthal asymmetry where the investigated 
effects, such as the elliptic flow coefficient V2, are of simi- 
lar order. For the case of the wounded- nucleon model [U 
and for the binary collisions the situation is illustrated in 
Fig. [1] The picture on the left shows all nucleons in both 
nuclei, the middle one the wounded nucleons, and the 
right one the binary collisions. We notice the mentioned 
effects for the distributions of the wounded nucleon and 
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the binary collisions: the twist of the principal axes, de- 
noted by the skewed lines, and the displacement of the 
center of mass, represented by a dot at the intersection 
of the principal axes. Statistical analyses may be carried 
out in the reference frame fixed by the reaction plane 
(we call it fixed-axes), or (in each event) in the frame 
defined by the shifted and twisted principal axes of the 
quadrupole moment (we call it variable- axes 1 ). 

The first part of the paper discusses the fixed-axes 
and variable-axes harmonic moments and radial pro- 
files obtained numerically from the Glauber Monte Carlo 
studies in several models: the conventional woundcd- 
nucleon model [Til ] , a model admixing binary collisions to 
wounded nucleons [12l. [l~3T] . a model with hot spots, as well 
as the hot-spot model where the deposition of energy oc- 
curs with a given probability distribution fSect. UTTjl . The 
results are presented in Sect. HVI andfVl The main result 
here is that the fixed-axes quadrupole moments, s, and 
their scaled standard deviation, Ae/e, vary significantly 
from model to model. The same holds to a lesser extent 
for the variable-axes moments, e*. On the other hand, 
the dependence of the scaled standard deviation Ae* /e* 
on the chosen Glauber-like model is weak, at most at the 
level of 10-15% for intermediate impact parameters. For 
all considered models the values range from about 0.5 for 
central collisions to about 0.3-0.4 for peripheral collisions. 
We examine the dependence on the mass number, provid- 
ing results for the gold-gold and copper-copper collisions. 
We also investigate the effects of the assumed weighting 
power of the transverse radius in the definition of the har- 
monic moments, finding that the choice is not important 
for studies of fluctuations. 

In Sect. I VII we examine the role of the center-of-mass 
and quadrupole-axes fluctuations on jet quenching. Ex- 
cept for very central collisions, the effect of the increased 
eccentricity of the opaque medium is canceled by the shift 
of its position and axes rotation, leading to almost no 
change in the azimuthal asymmetry of the jets leaving 
the interaction region. 

In Sect. I VIII we argue that the variable- axes quanti- 
ties are dominated by sheer statistics and certain prop- 
erties of variable-axes distributions can be explained 
in an elementary way through the use of the central 
limit theorem. In particular, in the absence of corre- 
lations between the location of sources and for central 
collisions we get the result of an appealing simplicity, 
namely Ae*/e*(b = 0) = y/i/n — 1 ~ 0.52, independent 
of the number of sources in the assumed model, the mass 
number of the colliding nuclei, or the collision energy. 
This result is fulfilled to a very good accuracy in actual 
numerical studies, where some correlations are present. 
For non-central collisions appropriate expansions are pro- 
vided. We also analyze the variable-axes profiles in this 



way. The effects of correlations between the location of 
sources are discussed in Appendix iDl 

In Sect. IVIIII we propose another method of encoding 
the information on the initial state, where each harmonic 
(including the odd ones) is evaluated in its own eigen- 
axes. The method can be used as a base for a smooth- 
ing procedure in preparation of the initial conditions for 
event-by-event hydrodynamic studies. 

In Sect. IIXI we make several comments referring to 
the collective flow. We note that the statistical anal- 
ysis of the variable-axes parameters e* carries over to 
the analysis of the variable-axes elliptic- flow coefficient, 
d|. For central collisions (in the absence of correla- 
tions) we find Av2/v%(b = 0) = ^4/tt- 1 ~ 0.52, inde- 
pendently of multiplicity, mass number, or the collision 
energy. This value is in the ball park of the recent ex- 
perimental data [1, Moreover, under the assump- 
tion of smoothness that most likely holds in hydrody- 
namics, which allows for perturbation theory around the 
azimuthally symmetric solution, one obtains the relation 
v 4 ~ v 2 2 f° r the octupole flow coefficient. Consequently, 
for the event-by-event fluctuations we find the prediction 
Avl/vl = 2Av* 2 /v* 2 . 

Appendices contain some more technical material, in- 
cluding the derivations of the statistical formulas. A sim- 
ple one-dimensional toy model illustrating the essence of 
the statistical intricacies is given in Appendix [Cj 



II. NOTATION 

In our study we use the standard Woods-Saxon nuclear 
density profile for the nucleus of mass number A, 



where the constant c, given in Appendix IB1 is such that 
the normalization J Airr 2 dr n(r) — A is fulfilled. For the 
considered gold and copper nuclei the parameters are 

# = 6.37 fin, a = 0.54fm, ( 197 Au), (2.2) 
i? = 4.14 fin, a = 0.57fm, ( 62 Cu). (2.3) 

A popular way to simulate the short-range repulsion in 
Glauber-like calculations 2 is to enforce that the centers 
of nucleons in each nucleus cannot be closer to each other 
than the expulsion distance of d = 0.4 fm. This feature 
is simple to implement in Monte Carlo generators. Some 
details are provided in Appendix [A"l 

We use the following standard convention for the axes 
of the reference frame: the z-axis is along the beam, 
the x-axis lies in the reaction plane, and the y-axis is 



1 We find this nomenclature more descriptive than the terms stan- 
dard and participant used in the literature. 



2 We note that this repulsion increases slightly the size of the nu- 
cleus, but the effect is negligible, see Appendix El 



FIG. 1: (Color online) Snapshot of a typical gold-gold collision in the x — y plane, 6 = 6 fm. Red and black circles indicate 
nucleons from nuclei A and B, respectively, plotted with the size (|3.2|l . The left picture shows all nucleons, the middle - the 
wounded nucleons only, and the right - the centers of mass of pairs of nucleons undergoing binary collisions. The straight lines 
indicate the (twisted) principal axis of the quadrupole moment, the blue dots show the center of mass of the system, while the 
outer circles denote the Woods-Saxon radius of gold, R = 6.37 fm. The units on the x and y axes are fm. 



perpendicular to the reaction plane. The azimuthal an- 
gle (j> S [— 7r, 7r] is measured relative to the y-axis, thus 
y = pcos(j), x = /?sin</i, where p is the transverse radius. 

We refer to the analysis in the fixed reference frame 
of the reaction plane as fixed-axes (sometimes called 
standard in the literature), and to the analysis where the 
particles in each event are translated to the center-of- 
mass frame and aligned with the major principal axis of 
the quadrupole moment as variable-axes (also called 
participant) . 



III. MODELS 

We describe briefly the models studied in this paper. 
The standard implementation of the wounded nucleon 
model at RHIC energies assumes that the inelastic cross 
section of the nucleon is 

a w = 42 mb. (3.1) 

The nucleon from one nucleus gets wounded when it 
passes closer to a nucleon from the other nucleus than 
the hard-sphere radius 

' ~~ (3-2) 



Then the weight w = 1/2 is attributed to the point in the 
transverse plane at the position of the wounded nucleon. 
The weight can be thought of as a measure proportional 
to the amount of the deposition of the transverse energy, 
which then is carried away by the produced particles. For 
studies of fluctuations only the relative weights are im- 
portant, and the overall normalization of the total weight 
can be chosen arbitrarily. In what follows we renormal- 
ize the distributions in all models to the number of the 
wounded nucleons, N w . 

For binary collisions the weight w = 1 is attributed to 
each collision point, which is taken as the mean of the 
coordinates of the two colliding nuclei. 



A successful description of multiplicities at RHIC 
has been achieved with a mixed model, amending the 
wounded nucleon model [ll[ with some binary colli- 
sions [HI, [l3|]. In this case a wounded nucleon obtains 
the weight w = (1 — a)/2, and a binary collision the 
weight w = a. The total weight averaged over events 
is then (1 — a)iV w /2 + aiVbin- The fits to particle mul- 
tiplicities of Ref. [l3[ give a — 0.145 for collisions at 
V^vw = 200 GeV, and a = 0.12 for ^fs^N = 19.6 GeV. 

Next, we consider a model with hot spots in the spirit 
of Ref. [HI , assuming that the cross section for a semi- 
hard binary collisions producing a hot-spot is small, 
"hot-spot = 0.5 mb, but when such a rare collision occurs 
it produces on the average a large amount of transverse 
energy equal to acr w /(7hot-spot- 

Each source from the previous models (wounded nu- 
cleon, mixed, or hot-spot) may deposit the transverse 
energy with a certain probability distribution. To incor- 
porate this effect, we superimpose the T distribution over 
the distribution of sources, multiplying the weights of the 
considered model with the randomly distributed number 
from the gamma distribution 



](w,k) 



k k exp(— kw) 



(3.3) 



This distribution gives the average value equal to w — 1 
and the variance var(w) = 1/k. This is at no loss of 
generality, since, as already mentioned, the individual 
weights used for carrying the statistical averages can be 
normalized arbitrarily. In this paper we do this superpo- 
sition on the hot-spot model, where the considered effects 
are largest. Thus, we take the weights (1 — a)g(w, k)/2 
for the wounded nucleons and ag(w, K)cr w /o'hot-spot for 
the binary collisions. We take k = 0.5, which gives 
var(iu) = 5. We label this model hot-spot+T. 

The four considered models, wounded nucleon, 
mixed, hot spot, and hot-spot+r, differ by the num- 
ber of sources and the amount of the built-in fluctuations. 
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For instance, the superposition of the T distribution with 
low values of k increases the variance. This increase is 
also generated with hot spots, which effectively reduce 
the number of sources. All these effects will be studied 
in detail below. 



IV. FIXED-AXES HARMONIC MOMENTS 

When the reaction plane is determined for each event 
(which of course can never be achieved exactly in the 
experiment, see e.g. Refs. [HI, [HI]), one can then choose 
the reference frame fixed by the reaction plane. The two- 
dimensional (boost-invariant) profile of the density of 
sources, /(/?, <j>), is obtained by averaging over the events 
belonging to a particular centrality or impact parame- 
ter class. The symmetry f(p, </>) = f(p, n — <j>) excludes 
odd components in the Fourier decomposition, while for 
equal colliding nuclei the symmetry f(p, <fi) = f(p, — 0) 
eliminates the sin(/(/>) functions. Thus, 

f(p, <t>) = fo(p) + 2/ 2 ( P ) cos(20) + 2f A (p) cos(#) + . . . , 

(4.1) 

where p is measured from the center of the geometric 
intersection of the two nuclei. The harmonic moments 
obtained form (|4.1[) . which we call fixed-axes, are also 
called "standard" in the literature. 

We first have a look at the Fourier profiles /; (p) , where 

/ = 0, 2, 4, In Appendix [Bl we show that at low values 

of p 

fl(p)~/f, (p<6), (4.2) 
while at high values of p 

Mp)~exp(-2p/a)p 2 (^j , (p»6). (4.3) 

Such a behavior is typical of Fourier expansions of 
smooth functions. 

We present the profiles fi(p) obtained with Monte 
Carlo simulations for the considered models in Fig. [2] 
The distributions in all models are normalized to the 
number of wounded nucleons, i.e. 

J pdpd^fip, <j>) = j 2irpdpf a (p) = N w . (4.4) 

We note that the profiles are substantially different from 
model to model. The monopole profile /o is broadest in 
the wounded nucleon model, then passing through the 
mixed model and the hot-spot model we arrive at the 
hot-spot+r model, which is most sharply peaked at the 
origin. Correspondingly, the quadrupole profiles f%{p) 
are concentrated further or closer to the origin. We note 
that the amplitude of subsequent harmonics decreases 
fast, such that taking I up to 4 is sufficient for any prac- 
tical matter. 



— wounded 

— - mixed 
hot-spot 
hot-spot + r 




hot-spot 
hot-spot + r 




hot-spot 
hot-spot + r 




FIG. 2: (Color online) The fixed-axes profiles fi(p) for gold- 
gold collisions in the analyzed models at several values of the 
impact parameter: b — 0, 4, 8 fm. The legend explains the 
assignment of line-types to the models. For 6 = 4 and 8 fm 
for each model there are three curves, corresponding to i = 0, 
2, and 4 from top to bottom. The i = 4 component is tiny. 
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In order to have some convenient quantitative mea- 
sures of the profiles of Eq. (|4.ip one introduces their ra- 
dial moments 



£fc,2 



f2npf l (p)p k dp 
J 2irpf (p)p k dp 



= (4-5) 



where we have introduced the moments 
1 

N w Jri 



h,i = -T7- / 2irpdpf l (p)p k 



(4.6) 



for a future reference. The choice of the weighting power 
k is arbitrary, with the typical choice k = 2. Higher val- 
ues of k make the measure more sensitive to the outer 
region of the system. We note that in the popular nota- 
tion 



£std — £2,2 



: s. 



(4.7) 



We observe that in all the considered Glauber models 
£ is practically independent of the model (top panel of 
Fig. [3]). Tiny differences come from different distribu- 
tions of the wounded nucleons and binary collisions. On 
the other hand, the scaled standard deviation, shown at 
the lower panel of Fig. [3J displays a strong dependence 
on the model at low values of b, with the hot-spot+r 
model yielding about twice as much as the mixed model. 
We also notice a strong dependence on b. At b = the 
curves diverge, which is an artefact of dividing by the 
vanishing value of e. The fluctuations are larger in mod- 
els effectively having the lower number of sources, which 
is obvious from the statistical point of view. 

As already noted in Refs. [T3, [Hj], the value of e ob- 
tained with the color glass condensate (CGC) is substan- 
tially higher than in all Glauber-like models analyzed in 
this paper. For comparison, the CGC result is shown as 
the upper curve in the top panel of Fig. [3j After the 
e-print version of this paper has been posted, a calcu- 
lation of fluctuations of e* in the CGC framework has 
appeared p^ |. The results are overlayed in in the bot- 
tom part Fig. [7j We note that at intermediate values of 
b the CGC values of Ae*/e* are significantly lower than 
in the considered Glauber models. 

In Fig. 2] we show the results for the octupole mo- 
ment, £4,2, and its standard deviation, obtained in the 
wounded-nucleon model. We note a very flat shape of 
£4.2 ~ b A at low b. This behavior is a direct consequence 
of integrating Eqs. (|4.2I4.3[) . The standard deviation 
grows with increasing impact parameter. In other mod- 
els considered in this paper the results are qualitatively 
similar. 



V. VARIABLE-AXES HARMONIC MOMENTS 

As is well known, the reaction plane cannot be de- 
termined precisely in an experiment, or not at all, 
which originated multiple methods of analyzing az- 
imuthal asymmetry in heavy- ion collisions. As has re- 
cently been realized P, 0, H, 0, S Hi , the purely statistical 
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FIG. 3: (Color online) The harmonic moment e = £2,2 and its 
scaled standard deviation for the analyzed models plotted as 
functions of the impact parameter. The result for the color- 
glass condensate comes from Ref. [lj. Gold-gold collisions. 



fluctuations caused by the finite number of particles lead 
to sizeable effects of the variable geometry in the initial 
stage of the collision. The effect can be seen qualita- 
tively in Fig. [I] where we notice a highly irregular shape 
of the distributions of both the wounded nucleons and 
the binary-collisions distribution. The mere presence of 
the fluctuation of the initial condition is obvious. What is 
somewhat surprising, however, is its size, leading to no- 
ticeable effects in the analysis of azimuthal asymmetry 
even at large numbers of participating nucleons. 

The center of mass of the distribution of the sources 
is not located at the geometrical center of the overlap of 
the two-colliding nuclei. In each event it is shifted in the 
x and y directions by Ax and Ay, defined as 



Ax = 



Ay 



(5.1) 
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FIG. 4: (Color online) The fixed-axes octupole moment, 
£4 = £2,4, and its standard deviation, plotted as functions 
of the impact parameter. Wounded nucleon model, gold-gold 
collisions. 



Ax, Ay [fm] 




FIG. 5: (Color online) The root mean square shifts of the 
center of mass in the in-plane direction, Ax (solid lines), and 
in the out-of-plane direction, Ay (dashed lines). The lower 
lines are for the wounded-nucleon, the middle for the hot-spot, 
and the upper for the hot-spot+r model. 

where i labels the sources and Wi denotes the weights. 
For an uncorrelated distribution of a large number of 
sources n one has 

(Ax) 2 = -R 2 x (w 2 ), (Ay) 2 ^-R 2 y (w 2 ), (5.2) 

where 

Rl = ^- J f(p, <j>)p 3 sin 2 <j)dpd(j), 

R l = J^J /(P. cos2 <M^' ( 5 - 3 ) 

are the mean squared radii of the geometric fixed-axes 
distribution, and (w 2 ) = J dw w 2 g(w, k). Since our sys- 
tem exhibits some correlation between the location of 




p[fm] 



FIG. 6: (Color online) The variable-axes profiles f*(p) for the 
analyzed models at several values of the impact parameter. 
The solid, dashed, and dotted lines correspond to i = 2, 4, 
and 6, respectively. The four models: hot-spot+r, hot-spot, 
mixed, and the wounded nucleon model yield curves arranged 
from the top to bottom, respectively. Gold-gold collisions. 
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sources, the formula (|5.2[) is not realized exactly, but hold 
qualitatively. In Fig [5] we show Ax and Ay as functions 
of the impact parameter for the wounded-nucleon (lower 
curves), the hot-spot (middle curves) and the hot-spot+r 
(top curves) models. The shift is more important for 
peripheral collisions, where n is lower, even though the 
source size itself decreases. 

The shift of the center of mass is physically relevant 
in the jet analysis, because it moves apart the formed 
fireball from the jet production points (Sect. rVTj) . 

In each event one can compute the principal axes of the 
ellipse of inertia. This corresponds to the choice k = 2 
as the weighting power. The angle between the major 
half-axis of the ellipse and the y axis is given by 



tan(20*) = 2 



(xy) - (x)(y) 
var(y) — var(x) 



(5.4) 



where the brackets denote the averaging over the parti- 
cles within the given event. Importantly, in the variable- 
axes calculations all coordinates in the given event are 
always shifted by (Ax, Ay) such that (x) = (y) = 0. 

The angle 4>* fluctuates sizably from event to event. In 
our notation, the superscript * indicates quantities aver- 
aged in such a way, that first in each event the rotation 
angle <ff is determined according to Eq. (|5.4[) . then the 
rotation is performed to the current principal-axis sys- 
tem, and finally the summation is done. As a result, 



/*(P,0) = /o(p) + 2/ 2 *(p)cos(20- 
+ 2/ 4 *(p)cos(40-40*)4 



(5.5) 



We call the above profiles the variable-axes profiles. Ob- 
viously, /o(p) = fo(p) and ft{p) > f t (p) for each p. 

In analogy to Eq. (|4.5[) we introduce the variable-axes 
moments 



-fcj 



J2npf l *(p)p k 
J2npf*(p)p^ 



(5.6) 



In the common notation for the variable-axes or partici- 
pant deformation parameter one has 



-part 



£ . 



(5.7) 



The profiles and moments for higher harmonics are sup- 
pressed, similarly to the fixed-axes case. This is clear, as 
the higher harmonics are evaluated relative to the axes 
determined by maximizing the quadrupole moment. As 
a result, only a few moments are needed to effectively 
parameterize the profile. 

Figure [H] shows the variable-axes profiles for I = 2,4,6 
(the I = profiles are equal to the fixed-axes case) for all 
considered models. The line-types distinguish the Fourier 
label I, while the hot-spot+r, hot-spot, mixed, and the 
wounded nucleon models yield curves for each I arranged 
from the top to bottom, respectively. Comparing Figs. [2] 
and [6] we note sizeable departures of the variable-axes 
profiles /* from the fixed-axes profiles fi (I = 2,4, ... ), 
in particular at small values of the impact parameter. For 
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FIG. 7: (Color online) The harmonic moment e* = €2,2 an d 
its scaled standard deviation for the analyzed models plotted 
as functions of the impact parameter. Gold-gold collisions. 
The open circles in the bottom figure show the calculation in 
the color-glass-condensate framework taken from Ref. [ljj. 



the central collisions (b = 0) the variable-axes profiles are 
non-zero solely due to fluctuations, as will be discussed 
in Sect. IVLT1 

In Fig. [7] we show the quadrupole moment e* and its 
scaled standard deviation. We observe a strong model 
dependence of e* at low values of b, with models having 
effectively lower number of sources yielding higher values. 
At b = the hot-spot+r model yields three times more 
than the wounded nucleon model. For all models the 
scaled standard deviation is close to the value 0.5 for 
central collisions and drops to about 0.3 at b = 14 fm. We 
argue in Sect. lVIll whv the central value is always close to 
0.5, independently of the effective number of sources. At 
intermediate values of b the relative difference in Ae* /e* 
between various considered models is at the level of 10- 
15%, which is not a very strong effect. 
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TABLE I: Dependence of the quadrupole asymmetry param- 
eters on the weighting power k from definitions (|4.5I5.6[1 



k 


2 4 6 


b = 


£ k,2 


0.047 0.064 0.089 0.121 
0.53 0.52 0.52 0.52 


b = 8 fm 


Ek,2 

Ae fc , 2 /£fc,2 


0.147 0.278 0.388 0.466 
0.44 0.48 0.48 0.49 


e k,2 

A £ fc,2/ £ fe,2 


0.176 0.319 0.452 0.555 
0.44 0.38 0.34 0.31 



We have examined numerically the role of the weight- 
ing power of the radius, k, entering the definitions 
(I4.5I5.6[) . This is important for the method, since as 
pointed out earlier, the value of k is arbitrary. The re- 
sult is that both £^2 and e£ 2 increase substantially with 
k, however the scaled standard deviation remains quite 
stable, in particular at low impact parameters. The re- 
sults are collected in Table [I] We vary k between and 
6, which is a wide range. At b = the scaled standard 
deviation remains practically constant, while at b — 8 it 
varies by 10% for the fixed-axes case and 25% for the 
variable-axes case. Due to this rather weak dependence, 
the particular choice of the weighting power k is not es- 
sential in studies of this quantity in event-by-event fluctu- 
ations. However, the values of £ and e* itself are sensitive 
to the choice of k. As already mentioned, the higher val- 
ues of k increase the sensitivity to the profiles at higher 
values of the transverse radius p. 

In Fig. [5] we present the variable-axes harmonic mo- 
ments and their scaled standard deviation for the copper- 
copper collisions. Due to the much lower number of 
sources compared to the gold-gold case of Fig. we note 
higher values of e* at low b. On the other hand, the scaled 
standard deviation is remarkably similar to the gold-gold 
case, especially at low b, as should be according to the 
arguments of Sect. lVIIl Certainly, the dependence on the 
mass number is a sensitive probe of the whole approach. 

In Sect. IvTlI we will show that many of the qualitative 
and quantitative features of the Fourier distributions as 
well as their moments have simple explanations on purely 
statistical grounds. 



VI. JET QUENCHING 

Jet quenching in dense matter occurs mainly at the 
very first stages of the collision [2(1 Hl| . The values of 
the nuclear modification factor 



Raa(pt) 



dN AA 
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measured in central Au+Au collisions for pt > 3 GeV 
fall significantly below 1. The dependence of the nuclear 
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FIG. 8: (Color online) Same as Fig. [7] for the copper-copper 
collisions. 



modification factor on centrality can be understood as 
due to the change of the size of the opaque medium which 
modifies the mean length of the path of the jet in the 
fireball. On the other hand the azimuthal asymmetry of 
the high-pT particles is believed to be a consequence of 
the geometric eccentricity of the medium. The difference 
of the path lengths for the jets moving "in plane" and 
"out of plane" leads to an asymmetry in the jet energy 
loss [H, [23|, [HJ . The angle-averaged nuclear modification 
factor Raa is very weakly dependent on the shape of the 
opaque medium, once its size has been fixed. 

In this paper we use a simple model 0, [2|| of the 
energy loss in order to explore the role of the shape of 
the event-by-event rotated absorbing medium. Neglect- 
ing the transverse expansion of the fireball for early times 
relevant for jet quenching, we expect that the shape of 
the medium is close to the initial conditions as obtained 
from the distribution of sources in the transverse plane 
from the Glauber models. Several prescriptions for the 
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distribution of the density in the transverse plane have 
been used, such as the wounded nucleon density, the bi- 
nary collisions density, the mixed density, or the color 
glass condensate estimate. All of these approaches use 
an event-averaged shape of the medium in which the jets 
propagate. However, it is clear that in each particular 
event the thermalized dense medium has a slightly dif- 
ferent shape and position with respect to the geometric 
reaction plane. In order to take this into account, one 
can use the variable-axes density f*(p, <p) as the density 
of the scattering centers for the propagating parton. The 
resulting increase of the eccentricity of the medium is 
expected to increase the asymmetry of the jet absorp- 
tion. A very similar effect has been discussed for source 
the profiles calculated including saturation in the CGC 
model. Drescher et al. [26J have found an increase in v 2 
by about 10 - 15%. 

The partons are produced in p-p collisions with the 
power-law spectrum dN/dpj, oc l/p^ 1 (the fragmenta- 
tion is not included), and the energy loss is taken as 

AE = fiE / M—^— f*(x +v x (l+T ),yo+Vy(l + To)), 
Jo 1 + T o 

(6.1) 

where the jet production point (xq, yo) is generated from 
a binary collision in the fixed-axes frame. The fact that 
we must use the fixed-axes frame here results from an 
absence of correlations of the very rare jet-production 
collisions and the soft collisions generating the opaque 
medium. In Eq. (|6.ip the initial time is denoted as to, 
while the time measured from tq is I. The direction of the 
parton transverse velocity (v x , v y ) is chosen randomly. 
For each choice of the model of the opaque medium the 
parameter /i in the energy loss formula (|6.ip is fitted to 
reproduce the high-p^ nuclear modification factor Raa 
|27| . Then the elliptic flow coefficient is calculated at 
different centralities (v 2 and Raa are pr-independent in 
such a simplified model). The variable- axes medium has 
a different shape from the fixed-axes one, with a larger 
eccentricity. As has been noticed, the raise in the geo- 
metrical eccentricity increases the asymmetry of the en- 
ergy loss [IB] . The variable-axes medium in the hot-spot 
model has an eccentricity of about 0.4, and the CGC 
calculation gives 0.5 at intermediate impact parameters, 
therefore one would expect a similar increase in v 2 at high 
Pt- However, there is one important effect that should 
be taken into account for the event-by-event modified ab- 
sorbing medium. The absorbing medium formed in each 
event is rotated and also shifted. The shift with respect 
to the fixed-axes frame is quite important (c/. Fig. [5]), 
yielding about 1/3 of the total effect. 

The resulting elliptic flow (Fig. [9]) at centralities larger 
than 20% resulting from the energy loss calculated with 
the wounded-nucleon model in the fixed-axes frame (solid 
line) comes out remarkably similar to the result of the 
hot-spot model in the variable-axes frame (dashed line). 
Only when the shift and rotation of the opaque medium 
are neglected (dotted line) the modification of the shape 
leads to an increase of the high pt elliptic flow coefficient 




FIG. 9: (Color online) Elliptic flow coefficient at high pr as a 
function of the number of wounded nucleons, obtained using 
the fixed-axes density of the wounded nucleons f(x, y) in the 
energy loss formula 16 . 1 1 (solid line), and with the variable-axes 
density f*(x,y) for the hot-spot scenario (dashed line). The 
dotted line represents the result for the variable-axes density 
but without the shift and rotation of the opaque medium. 

v 2 by about 10 — 15%, We have checked that the cancel- 
lation of the effects of the increased eccentricity of the 
medium and of the shift and rotation with respect to the 
jet emission points at larger centralities happens also for 
other considered models. 

In the forthcoming experiments at the Large Hadron 
Collider (LHC) the analysis of jet tomography with re- 
spect the the event- by- event reconstructed plane must 
take into account the relative shift and rotation effects 
discussed above. 



VII. STATISTICAL INTERPRETATION OF 
THE VARIABLE-AXES MOMENTS 

In this section we analyze the variable-axes moments 
and profiles from the viewpoint of statistical methods. 
The purpose of this study is to understand certain fea- 
tures of the numerical results presented earlier on more 
general formal grounds. It turns out that quite simple 
expressions can be found for the case where correlations 
between the location of sources are neglected, which al- 
lows for the standard usage of the central limit theorem. 
The Glauber models do induce some correlations, as can 
be seen from Fig. [TJ For instance, a nucleon from the 
skin of one of the nuclei, as present in the middle picture, 
wounds several nucleons from the other nucleus. As a re- 
sult, correlation between the locations of the wounded 
nucleons is generated. For simplicity, we neglect all such 
correlations in the analytic analysis of this section. Their 
role is discussed in Appendix [Dl 

If such correlations are strong, their analytic inclusion 
is difficult and one has to resort to numerical simulations 
such as those presented in the earlier sections. We also 
take all weights equal to unity, «;,• = 1, in order to avoid 
notational complications. 
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Appendix [C] contains a very simple analysis in a 
one-dimensional independent-particle toy model, which 
avoids some notational complications but grasps all es- 
sential features of the full case. 

From definition of the variable- axes moment, we need 
to evaluate 



~-k.i 



3 = 1 l J 3 



n 

— «-£p£cOS[J(&-^)]», 

4,0 n ^ J 



(7.1) 



3=1 



with ((.)) denoting the averaging over (infinitely many) 
events and j labeling the source. The weighting, accord- 
ing to the definition (|4.5p . is done with the transverse 



radius to the power k, i.e. 
notation 



Y, 



1 



3 = 1 



p$ cas(l(/>j), Xi 



Let us introduce the 



^sin^) (7.2) 

3 = 1 



(recall that we measure the azimuthal angle from the y- 
axis). The rotation angle (jf depends on the distribution 
of particles in the given event, by definition maximiz- 
ing the quadrupole moment (I — 2), i.e. the quantity 
- Y^j=x Pj cos[2(0j —(/)*)]■ This gives the relations 



cos(2(f ) 
sin(2cf ) 



(7.3) 



Let us denote Xj — (pj , <f>j ) as the short-hand nota- 
tion for the polar coordinates of the source point, and 
F(xi, x%, . . . , x n ) as the n-particle probability distribu- 
tion. Then we can rewrite Eq. (|7.ip as 



- k.l 



— !— [ dxi . . . dx n F(x 1 , . . . , x n ) x (7.4) 

Ikfi J 
1 - 

- pj[cos(l4>j) cos(2(f)*) - sm(%) sin(20*)] = 

Y t Y 2 + X t X 2 



3 = 1 



1 f 

- — / dxi . . . dx n F(xi, x n )- 

J-kfi J 



Let us analyze the quadrupole moment e* k 2 , which is 
the simplest but also the most important measure. For 
that case Eq. (|7.ip becomes 



'fc,2 



- — [ dxi . .. dx n F(x ll 

J-kfi J 



Yi+Xl 



— (( X /Yi+XD) 

1 kfi 



(7.5) 



nI k .o 



« 



\ 



J2 Pj cos(20j 



£p*siii(2&) | )). 



Thus, the variable-axes quadrupole moment corresponds 
to a highly "non-local" average, involving infinitely many 
moments through the square root function. 

In the absence of correlations between collision points 
the many-particle probability distribution factorizes into 
F(xi, . . . ,x n ) = f(xi) . . . f(x n ), where f(xi) are normal- 
ized to unity. Hence all the information on the system is 
contained in the one-particle distribution functions (|4.1[) . 
or, equivalently, in the fixed-axes profiles fi{p). The 
variable-axes profiles /;*(p) and moments e* kl are then 
expressible in terms of the fixed-axes quantities. In or- 
der to get some useful relations, we use a method similar 
to the techniques of Refs. 0, H SI]. In the limit of 
large n the goal is accomplished with the help of the cen- 
tral limit theorem. Indeed, the variables p k cos(2</>) and 
p k sin(2</>) entering Eq. (|7.5p are independent from one 

another, since / p 2k cos(2^>) sin(2</>) = 0. Therefore for 
sufficiently large values of n one may use the central limit 
theorem, implying the normal distribution for the vari- 
ables Y2 and X2 ■ The details of this calculation are given 
in Appendix [Pi The final results can be cast in the form 
of a series involving the confluent hypergeometric func- 
tions: 



_ ya V (26a 



2 \m 
Y 2 ) 



(7.6) 



r (m + |) T (m + § ) i*i (-±;m + 1; -Jh) 



|2 



where the average and standard deviation of the variables 
(|7.2j) are (see Appendix iDjl 



(7.7) 



% = 


Ik,2, 


< = 


—-{hkfi - 27fe 2 + hk.A) 




—-{hkfi — hk,i)t 




1 1 


6 = 


24 2<4 2 • 



For the special case of central collisions (where 6 = 0) 
only the m = piece contributes to the series (|7.6[) and 
we have the very simple result 



-k,l 



2Ik,oVn 



(6 = 0). 



(7.8) 



Similarly, for the scaled standard deviation in central col- 
lisions we obtain (see Appendix [D|l 



Ae 



kj 



-k.l 



- 1 ~ 0.523, 



(6 = 0) 



(7.9) 



Although the above result is approximate, as it has been 
obtained with the assumption of no two-particle corre- 
lations between the location of sources, it shows an im- 
portant feature present in the Glauber simulations. The 
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value of the scaled standard deviation in central collisions 
is close to 0.5 and is asymptotically independent of n (cf. 
Figs. [7] and [8]). The behavior is also seen in the numbers 
given in Table UJ for 6 = 0. 

One may also derive the expression for the profile func- 
tion f^ip) hi the large- n limit and for the case with no 
correlations. For central collisions the result is (see Ap- 
pendix |D]) 

r 2 { P )^\ /Z-_>/o(p), (6 = 0) (7.10) 

Thus the variable-axes quadrupole profile is proportional 
to the monopole profile times p k . Note that this behav- 
ior reflects the chosen weighting power fc, hence in this 
sense is technical rather than physical. One verifies that 
Eq. (|7.10p reproduces immediately Eq. (|7.8p . 



VIII. MULTIPLE-AXES PROFILES 

In Eq. (|5.5p the rotation angle </>* has been fixed with 
the quadrupole moment. All higher harmonics were ob- 
tained with respect to this angle. One can give another 
prescription, which is superior for encoding the shape of 
the system for studies of event-by-event fluctuations. For 
a distribution of sources in each event we may provide the 
harmonic profile, as well as its orientation relative to the 
y axis. Thus, each profile /* has its own rotation angle 
4>f ■ We introduce the expansion (the meaning of * is now 
different than in Sect. fVl and IIX[) 

f*(p,4>) = fo(p) +2/ a *(p)co8(2^-2&) (8.1) 
+2f*(p) cos(30 - 3$) + 2/ 4 *(p) cos(40 - 4#) + . . . 

This expansion is complete, as 

cos(l(j) — l(f>i) — cos(l(f>) cos(l(j)i) + sin(Z0) sin(Z</> ; *) = 
= ai cos{l<p) + bi sin(Z0), (8.2) 

which provides the full Fourier expansion of the distri- 
bution in each event, including both the sine and co- 
sine functions. Note also the presence of odd moments, 
I = 3, 5, ... . These moments average out to zero in the 
limit of infinitely many events, but have non-zero fluctu- 
ations from event to event. 

In Fig. [10] we show the multiple-axes profiles for the 
wounded nucleon model for the first few harmonics. We 
note that there is no longer a strong suppression as I 
is increased. Clearly, to describe completely the full 
shape of the distribution we need as many harmonics as 
sources! In actual applications certain smoothing must 
be included, which effectively cuts off the higher harmon- 
ics from the expansion. These equilibration processes 
cause sharp shapes to smooth out, i.e. high harmonics 
are damped. One may therefore propose a smoothing 
prescription based on multiple-axes profiles, introducing 
a suitable cut-off function in the Fourier index I. Details 




p[fm] 



FIG. 10: (Color online) The multiple-axes profiles f*(p), i = 
2, 3, 4, 5, 6, for the wounded nucleon model at several values 
of the impact parameter b. Gold-gold collisions. 
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FIG. 11: (Color online) The multiple-axes octupole moment, 
e% = £2,4, and its scaled standard deviation, plotted as func- 
tions of the impact parameter. Wounded nucleon model, gold- 
gold collisions. 



and application to event-by-event hydrodynamics will be 
presented elsewhere. 

In Fig. [11] we show the results for the multiple-axes oc- 
tupole moment, e\ 2 , and its scaled standard deviation for 
the case of gold-gold collisions. As explained above, this 
moment is computed in each event in the reference frame 
which maximizes the octupole distribution. We note that 
again at b = the scaled standard deviation is close to 
0.5, as in the large-n limit and in the absence of correla- 
tions it assumes the value As* k /e* k (b = 0) = ^/A/n — 1. 
The details of the analysis are given in Appendix [Dl 



IX. FLUCTUATIONS OF THE ELLIPTIC FLOW 

The statistical analysis for the fluctuations of the 
variable-axes shape parameters, in particular e*, carries 
over to the fluctuations of the elliptic flow coefficient v\* . 
These fluctuations, which are an important probe of the 
nature of the early-stage dynamics of the system [30l |. 
have recently been measured at RHIC [1, 0, 101. In fact, 
the experimental procedure used in these analyses iden- 
tifies the elliptic flow coefficient with the participant or 
variable axes V2, here denoted as u|- 

The relevance of studies of fluctuations of the initial 
shape of the fireball comes from the well-known fact that 
for small elliptic asymmetry one expects on hydrody- 
namic grounds the relation 



AvX 



Ae* 
1*~ 



(9.1) 



As argued in Ref. [31|, the result (|9.ip indicates that 
the mean free path in the matter created in the initial 
stages of the heavy-ion collisions is very small, although 
turbulence does not develop. 



The statistical method used in Sect. [V] carries over to 
i>2. An immediate consequence of Eq. (|9.ip , under the 
assumption of the absence of correlations between the 
location of sources, is the result for the variable-axes co- 
efficient, v%, in central collisions: 



(b = 0) 



- 1 ~ 0.52. 



(9.2) 



The values obtained in Refs. 0, [T(| for the scaled stan- 
dard deviation of v% are between 0.35 and 0.5 for all im- 
pact parameters. 

An argumentation for the result (|9.ip may be done 
on general grounds as follows: schematically one may 
denote the hydrodynamic equations as L(ip) = 0, where 
L is the operator for hydrodynamics (involving partial 
differentiation, etc.), and if) is the set of hydrodynamic 
functions of space-time describing the state of the system. 
If the evolution is smooth, one may expand to first order 
around the azimuthally-symmetric system tpo, 

L(V0 = L(^ + 6i>) a L(iM + L'iiMW, (9.3) 

where Sip is the asymmetric piece, and the prime de- 
notes the differentiation with respect to the hydrody- 
namic variables ip. Since L(V'o) = 0, we have to first 
order L'(tp )8ip = 0. Then, due to linearity of the equa- 
tion for 5ip, we have ~ II^VKMIL *- e -; the magni- 
tude of the solution at time t is proportional to the initial 
condition at to. This concerns all hydrodynamic proper- 
ties, in particular the shape and flow. As a result, the 

coefficient determined from the momentum spectra at 
time t is proportional to the initial spatial quadrupole 
asymmetry e*. The feature holds event-by-event, hence 
the result (|9.ip follows. Thus Eq. (|9.ip is a consequence 
of applicability of perturbation theory for the small de- 
parture from cylindrical symmetry. 

We note that Eq. (|9.ip holds separately for the fixed- 
axes and the variable-axes analyses, with the obvious re- 
quirement to use the same method on both sides of the 
equation. 

One may ask if a similar argumentation can be used 
for the higher harmonic flow coefficients, u|, etc. The 
results of the previous sections show a strong suppression 
of subsequent harmonic moments of the distribution of 
sources in the case of the fixed-axes and variable-axes 
analyses. This suggests the hierarchy 



ijj = i/j + XS^2 + X 6i/j 4 + ■ ■ 



(9.4) 



where A, typically of the order of a few percent, is 
the small expansion parameter, while the subscripts 
0, 2, 4, . . . label the harmonics. Expansion of the hydro- 
dynamic evolution to second order in A, again under the 
assumption of smoothness, yields now 



L(i/j) = L(^o)+Ai'(^ )^2 

+ A 2 [L'(Vo)^4 + L"(Vo) W 2 ) 2 /2] 



(9.5) 
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We carry out the perturbation theory extracting the 
second-order equation, 

L'(V>o)^4 = -L"(^ )(^ 2 ) 2 /2, (9.6) 

and note that the evolution of SipA is coupled to (81P2) 2 , 
which acts as a source term in the linear inhomogeneous 
equation for 8^4. 

Various harmonic components evolve hydrodynami- 
cally with different time scales. Let us denote r 2 as the 
characteristic time for the operator L'(ipo), or S?p2, and 
T4 as the characteristic time for the operator L"(ip ). If 

t 2 >t 4 , (9.7) 

then the time scale for the source term in Eq. (|9.6p is 
much larger than for the operator L"(ip ). In that case 
for t > t 

\\Mt)\\~\\SMt)\\ 2 ~\\$Mto)\\ 2 - (9.8) 

In words, at late times the octupole deformation is pro- 
portional to the square of the initial quadrupole defor- 
mation, and looses memory of the initial octupole defor- 
mation ||^4(£o)||- In particular, this means that 

v * 4 „ e* 2 ~ v* 2 2 . (9.9) 

In Ref. [HJ the variable v^jv\ has been suggested as a 
sensitive probe of the hydrodynamic evolution. Also, the 
simulations of Refs. [32|, [H| show that with increasing 
time the value of vi saturates (suggesting very large T2), 
while V4 quickly assumes the value proportional to w|, 
supporting the assumption r 2 S> T4 used in the above 
argumentation. The data of Refs. [34[ comply to the 
result (|9.9p . perhaps except for very low values of the 
transverse momenta. 

For the fluctuations one gets immediately from 
Eq. (£2) 

Avt AvS Ae* 

^=2—^ = 2—. 9.10 

v\ v% e* 

Relation (|9 . 1 0|) , if verified experimentally, would support 
the scenario of smooth hydrodynamic evolution with the 
mentioned hierarchy of scales. 

At sufficiently late times all deformations are deter- 
mined by the initial e*. This results in other relations, for 
instance for the azimuthal Hanbury-Brown-Twiss (HBT) 
correlation radius, -Rhbt (</>), one expects 

Ri ~ Rl, (9.11) 

where -Rhbt(^) = Rq + 2i? 2 cos(20) + 2R 4 cos(40) + . . . 

X. CONCLUSION 

We have presented a comprehensive study of the shape 
fluctuations in a variety of Glauber-like models. Here is 
the list of our main points: 



1. We compare four Glauber-like models, with dif- 
ferent degree of fluctuation: the wounded-nucleon 
model, the mixed model, the hot-spot model, and 
the hot-spot model with the superimposed T dis- 
tribution. 

2. We obtain numerically the fixed- axes and variable- 
axes harmonic profiles and analyze their moments. 
The variable-axes moments £*, and the fixed-axes 
scaled standard deviation Ae/e are sensitive to the 
choice of the variant of the Glauber model, while 
the Ae*/e* is not, changing at most by 10-15%. At 
intermediate values of b the results of the Glauber- 
like models for Ae*/e* lie significantly above the 
color-glass-condensate predictions of Ref. . 

3. We present expansions for the variable-axes mo- 
ments and profiles. These analytic formulas explain 
the features of the simulations, in particular, they 
show that at b = the multiple-axes scaled vari- 
ances are close to the value 0.5, insensitive of the 
model used, the mass number of the colliding nu- 
clei, or the collision energy. In essence, the behavior 
of the scaled variance, used as a popular measure 
of the event-by-event fluctuations, is governed by 
the statistics. 

4. Unlike the results of Ref. [26| which finds an in- 
crease of the jet elliptic flow u 2 at the level of 10%, 
we find that the effect of the increased quadrupole 
eccentricity is largely canceled by the shift of the 
center of mass and rotation of the axes of the 
absorbing medium. This leads to practically no 
change of the jet emission asymmetry at inter- 
mediate and large impact parameters. Only for 
small impact parameters the appearance of the 
quadrupole moment in the shape of the medium 
wins over the relatively less important shift and ro- 
tation. 

5. We propose to use an improved harmonic expan- 
sion, the multiple-axes expansion, where the har- 
monics in each event are evaluated with their own 
reference frame. Such a scheme may be a starting 
point for the event-by-event hydrodynamic studies. 
The details will be presented elsewhere. 

6. The analysis of the variable-axes moments in the 
coordinate space directly carries over to the collec- 
tive flow and analysis of v\ in the momentum space. 
In particular, Eq. (|9.2[) holds for the variable-axes 
elliptic flow coefficient. 

7. Finally, we comment that under plausible assump- 
tions of smoothness, the hydrodynamic evolution 
leads to sensitivity of higher flow harmonics, U4, etc. 
to the initial quadrupole deformation only. Higher 
harmonics of the deformation are irrelevant. Then 
Eq. dHJ) or flgHUD follow. 
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APPENDIX A: DETAILS OF THE 
MONTE-CARLO PROCEDURE 

The 3-dimensional positions of the nucleons in a nu- 
cleus are randomly generated from the Woods-Saxon dis- 
tribution (|2.ip . Whenever the center of a nucleon is gen- 
erated closer than the expulsion distance of d = 0.4 fm 
to a center of any of the prior generated nucleons, this 
nucleon is discarded and generated anew. The procedure 
results in a certain "swelling" phenomenon. However, for 
the chosen value of d = 0.4 fm the effect is tiny, increasing 
the R parameter by 0.01 fm only, which is a small fraction 
of a percent. The swelling effect could be compensated 
by reducing appropriately the original R parameter of 
the distribution (|2.1j) . At higher values of the expulsion 
distance d the swelling effect increases. For instance, if 
one chose d = 1 fm, then the resulting size parameter 
is R = 6.575 fm, 3% larger than needed, and then the 
compensation in the original value of R should be done. 



APPENDIX B: PROPERTIES OF THE 
FIXED-AXES HARMONIC PROFILES 

At fixed values of b the shape of the profiles fi(p) 
reflects in a simple manner the average distribution of 
sources in a given model. Let us assume that d = 0, i.e. 
we ignore the short-range correlations, which complicate 
the analysis. The nucleus thickness function is 



T(s) 



dzn 



(Bl) 



with the density function n defined in Eq. (J2JJ) . The nor- 
malization is J 2irsdsT(s) = A. In the wounded nucleon 
model one has the following formula for the density of 
sources in the collision of large nuclei A an B: 

n w (b,p) = T A (p + b/2)[l-exp{-a w T B (p-b/2))] 
+ T B (p- 6/2) [1 - exp(-a w T A (p + 6/2))] . 

(B2) 

The corresponding formula for the binary collisions is 

7V bin (6,p) = a hin T A (p + b/2)T B (p-b/2). (B3) 

Let us introduce the short-hand notation u = p 2 + b 2 /4, 
v = pb sin <fi. We may then rewrite Eq. (|B2IB3|) in a more 



N w (u,v) = T A (u + v)[l — exp(—a w T B (u-v))] 
+ T B (u -v)[l — exp(—a w T A (u + v))] , 
N hin (u,v) = a hin T A (u + v)T B (u~v). (B4) 

When p « ii then v <C u and the low-p expansion of 
the source density corresponding to Eqs. (|B4|) has the 

form J2 n c n p 2n sin 2 ™ tf>. Because d^sin 2 ™ </>cos(2m</>), 
needed for the the decomposition (|4.ip . vanishes at 
n > m, we obtain the result (I4.2[) . 

The normalization constant in the Woods-Saxon func- 
tion (|2~T1) is 



-A/ 



.,(- 



,R/a 



(B5) 



Here Li„(z) = YlkLi z k /k n denotes the polylogarithm 
function. At high values of p the function (|2.ip asymp- 
totes to cexp ((R — r)/a)). Correspondingly, the nucleus 
thickness function at large s becomes 



T(s) - cV2 



iras e 



s/a 



(B6) 



For p 3> b we also have v <C u, hence with a calculation 
similar as for low p we find Eq. 



APPENDIX C: THE TOY PROBLEM 

This Appendix contains a detailed description of a 
toy model illustrating in a simple manner the statistical 
techniques used in the full-fledged calculation. Consider 
the one-dimensional problem (in the azimuthal angle (j)) 
where we randomly generate uncorrelated particles from 
a distribution containing the monopole and quadrupole 
moments only, 



r 1 X i 

£e[ -2<2 ] 



f{<j>) = l + 2ecos(2<£), 

The distribution has two fixed-axes moments, 
1 



(CI) 



/o = 



2tt 
1 

27 



#/(<£) = 1, 



#cos(20)/(0) 



(C2) 



Suppose we generate randomly n particles according to 
the distribution (|C1|) in each event, and subsequently 
carry the averaging of the results over the events, de- 
noted as ((.)). For instance, fi is estimated as 



1 ™ 

«-£cos(2&)», 



(C3) 



k=l 



where k labels the particles in each event. The equal- 
ity becomes strict as the number of events approaches 
infinity, which we assume implicitly from now on. 
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For the variable-axes quadrupole moment we need to 
rotate the particles with the rotation angle <jf , which 
changes from event to event. Thus we need to evaluate 
(f 2 has the meaning of e* from the other parts of the 
paper) 



1 n 



cosf2( 



')]»■ 



(C4) 



fe=i 



The rotation angle <jf depends itself on the distribution of 
particles in the given event. It is such that the quantity 
A = i cos [2(</ , fe — 4>*)] assumes maximum, which 

gives the conditions dA/d<j>* = 0, d 2 A/d(4>*) 2 < 0. The 
solution is 



cos(2<£*) = Y 2 /JY 2 +X. 



sin(2<n = X 2 /^Y 2 + X 2 , (C5) 
where we have introduced the short-hand notation 

^ n \ n 

Y 2 = - J2 cos(20 fe ), X 2 = - V sm(2</> k 



k=l fc=l 

Using the above formulas in Eq. (|C4I) yields 



/a = (i\/Yi+Xi)) 
= « 



(C6) 



(C7) 



We see that the variable-axes moment corresponds to an 
average of the square root of sums (|C5j) , thus is a highly 
"non-local" object. 

For sufficiently large multiplicity of the events, n, one 
may evaluate Eq. (|C7[) with the help of the central limit 
theorem. Consider the variables Ck — cos(2^>fe) and 
Sfc = sin(20fe). The average is 

Y n y r 27T 

5 = « - E c *» = ^ y = e - ( cg ) 



fc=i 



while for the variance we have 



(C9) 



1 

27 



fe=i 

2tt 



>cos(2^) 2 /(0)- e 2 = -- £ 2 . 



Likewise, for the Sk variable 



5 = «~E S *» = T- / #sin(20)/(^)=O, 
n fc=1 ^ j 

n , 



Importantly, there is no correlation between Y% and X 2 , 
as 

«- E = 5~ / #cos(20) sin(2^)/(</)) = 0. 



fe=i 



(Cll) 



According to the central limit theorem, the distribution 
of the Y 2 and X 2 variables is Gaussian, 



f(Y 2 ,X 2 ) = 



2ira c a s 
exp 



exp 



(Y 2 -c) 2 , X 2 



2ol 



2a\ 



ttVI - 2e 2 
Introducing the notation 



(Y 2 -e) 
l-2e 2 



■Xt 



(C12) 



Y 2 = q cos a, X 2 — qsina, q 2 — Y 2 + X 2 , 



S = 



1 



1 



1 



2a 2 . 2o 2 1 - 2e 2 



- 1, 



(C13) 



we may write 
exp 



ttVI - 2e 2 
q 2 + e 2 — 2gecosa 
l-2e 2 



(C14) 



nSq 2 sin 2 a 



What we will need below is the integral of this distribu- 
tion over a, 



daf(q,a) 



2n 



Wl - 2e 2 



exp 



E (« T{]+ . 1) ^ ( 2neq 



3=0 



J- 



l-2e 2 



As a check, it follows that 
qdqdaf(q, a) 



q 2 + e 2 
l-2e 2 

(C15) 
(C16) 



VT^2e 



E( 2e2 



3=0 



where we have use the definitions (|C13[) and the formula 



3=0 



31 



I- A 



(C17) 



We may now evaluate the variable-axes moment (|C4I) 
We have the following series: 



/a 



1 - 2e 2 °° 
q dq da qf(q, a) = — t=- V (2c 



2\3 



3=0 



r0- + i)r(i + |) 



, 1 ne 2 
lfl |-2' J + 1; "l^ 



fe=i 



(C18) 
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n=10, 20, 50, 100, 400 




FIG. 12: (Color online) Toy model. Dependence of the 
variable-axes moment /J on the fixed-axes quadmpole mo- 
ment e for several values of the number of particles n. As 
n increases, we pass from to to bottom with the presented 
curves. The straight line is the n —* oo limit, i.e. /| = e. We 
note that the effect of the departure of /| from e is strongest 
at low e and low n. 



f * 

E 7 



n=100 




FIG. 13: (Color online) Toy model. The rate of convergence 
of the series (|CT8)| . The curves from bottom to top we show, 
correspondingly, the results of summing up 1, 2, 3, 4, and 5 
terms, as well as the full result. 



We were not able to sum up this series into a closed form, 
but one may readily use it for practical calculations in a 
truncated form. At e = we have 



f*2 (e = 0) 



(C19) 



which shows the expected 1/ 1 ^fn behavior for a statistical 
fluctuation. 

The numerical results of the series (|C18[) are presented 
in Fig. [T2J We note that the effect of the departure of f% 
from e is strongest at low e and low n. 

Figure [13] shows the rate of convergence of the series 
(|C18p for n — 100, where we show the subsequent results 
of summing up 1, 2, 3, 4, and 5 terms. We note that 5 
terms are sufficient to achieve accuracy better than 1%. 

We note that an expansion of the result in powers of 
e is useless due to slow convergence properties. The first 
few terms are 



/a = 



1 



n — 1 



n 2 - Qn + 3 
18 ' 



O fn 3 e 6X 



(C20) 



CT ( f 2)/ f 2 n=10, 20, 50, 100, 400 
. 5 : 




FIG. 14: (Color online) Toy model. The dependence of the 
scaled standard deviation on e for several values of the number 
of particles n. The curves from top to bottom correspond to 
n = 10, 20, 50, 100, and 400, respectively. 

hence the effective expansion parameter is ne 2 . Thus for 
large values of n the convergence radius in e is very small. 

The evaluation of the second moment in the q variable 
yields 



VI - 2e 2 



{2e 2 ) J {(-2j + n-2)e 2 +j + l)r(j + ±) 



((Y 2 2 + X 2 )) = / qdqdaq 2 f(q,a) = 



E 

3=0 



1 + (n- l)e 2 



(C21) 



The obtained result is obvious from a direct evaluation 
form the definition. We compute 

((Y 2 + X 2 ))= (C22) 
<<(^f>s(2^ +fi|>n(2<fe)j )) = 

«~ + ~^ (c ^ 2 ^) cos(2^ fc ) + sin(2^-) sin(2<fc)))). 

Since the system is uncorrelated and (cos(20;)) = e, 
(sin(2</>/)) = 0, we immediately obtain 



((Y 2 +Xi)) = - + 



1 n(n — l)e 2 



(C23) 



n w 
in agreement with Eq. (|C21[) . 

From Eqs. (|C18IC21[) we may obtain the expression 
for the variance of the distribution of the variable-axes 
moment. A simple formula follows for the case e = 0, 
where 

2 



0.215 



var(/|) = i - ( £L) = LA^—. (C24) 
n \2y/n J n n 

The scaled variance and scaled standard deviation are 



var(/ 2 * 

fS 



2 \pa 

- V 0.242 



(C25) 



17 



Note that (for e = 0) there is no dependence on n in 
the scaled standard deviation. The general case of the 
dependence of the scaled standard deviation of / 2 on e 
for various values of n is shown in Fig. 1141 The result is 
obtained numerically from Eq. (|C18IC21|) . 

APPENDIX D: CENTRAL LIMIT THEOREM 
AND THE MULTIPLE-AXES MOMENTS AND 
PROFILES 

This Appendix contains some details of the applica- 
tion of the central limit theorem to the analysis of the 
multiple-axes moments and profiles. The calculation is 
carried out for the case where each harmonic moment has 



its own rotation angle, as described in Sec. IVIIIl Define 

n i n 

y ' = -E^ C0S (^)' x 1 = -J2p"Mi^) (Dl) 

" i=i U j'=i 

The rotation angle <fi* satisfies 



cos(Z0*) = Yl/^Yi 2 + Xf, 

sin(Z0*) = Xt/^Yf+X?, (D2) 

and it depends on the polarity index I. We need the 
averages 



Y, = 



X, = 



4 



'X, 



j dxi . . . dx n - ^2 Pj C0S ( l( Pj)f(xi, x n ) = / dcj)pdpf(p, <j>)p k cos(Z^) = J fc ,, 
J n _. =1 - j 

j dxi . . .dxn-'^p^ $m(l(j)j)f(xi, . . .,x n ) = / d<ppdpf(p,cf))p k sin(^) = 0, 
J n J= i J 

/ dxi . . .dx n — ^2 coa(l4>j) ^ p k , cos(l^)f(xi, . . .,x n ) - (Y{f 

■> n .7 = 1 .7'=1 



d(j)pdpp k cos(Z0)/(p, 4>)\ -q. 



= - j d^p 2fc cos(^) 2 /(p,0) + — 
n J n 

= kij d ^ d PP 2k f^ M 1 + «»(2ty)) - h 2 k , = i-(Z 2fe , - 2/2 , + J 2fe , 2i ) 

- OC OO 

/ efei . . . darn-g p* sin(^j) p£ sin(Z^/)/(a; 1 , . . . ,x n ) - (Xt) 2 
J n j=i j'=i 

If 1 

= - / d(f>pdpp 2k sin(7</>) 2 /(p, 4>) = — (7 2 fc,o " hk,2i)- 
n In 



(D3) 



r 



From the central limit theorem, the distribution of F and 
Xi has the normal form 



f{Y h X l ) = 



1 



2ir(TY l crx, 



exp 



(Xi-Yif X? 



24 



Introducing the q and a variables through 

Yi = gcosa, X; = qsina, 
we can rewrite Eq. (|D4j) as 

/(«,<*) 



exp 



2it(Jy 1 o'x 1 
Q 2 + q 2 , qq y cos a 



(D4) 



(D5) 



(D6) 



24 



+ <5g sin a 



'Y, 



where 



5 = 



2cr Y 2g\ 



Mill - hk,2l) 



(D7) 



(hk,0 — 2/ fe 1 + l2k,2l){hkfi — hk,2l) 



Next, we expand in the Taylor series in S and carry the 
integration over a, which yields 0, [l5j 



p2tt 

f(q) = da exp 
Jo 



^ (^ 2 sin 2 a) m 

7 r = ex P 

m=0 *' A! 

,M^/4)r(m + i) 



g 2 + (F;) 2 gF/ cos a 



2a 2 



g 2 + F; 2 



2a 2 



2£c4g 



m=0 



F 



(D8) 
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where I m is the modified Bessel function. The relevant 
moments of f(q) are 



qdqf(q) = 1, 



(D9) 



qdqqf(q) 



\/2er? °° 



'Y, \ - 



l-KOx. 



r (m 



r (m 



' m=0 



1*1 



(2^4 ) m x 



1;-t 



qdqq 2 f(q) 



hk,o + (n - 1)7? , 



where iiq is the confluent hypergeometric function. 

From Eq. (|D9[) one derives immediately the formulas 
for ^ = 2 listed in the main text as Eq. (|7.6l7.8p . For a 
general value of I the result for the variable-axes moments 
(recall each moment has its own rotation angle) reads 



£ k,l — 



(D1Q) 



r (m + i) T (m + |) 1F1 ( -|;m + 1; 



2ar 



Var ( £ fe,i) 



/2fe,0 + ("- 1 ) / fci 



nT 2 
'"fe,0 



(4«) s 



Figure [15] compares the formulas (|D10|) for the 
quadrupole case (/ = 2) with the Monte Carlo simulation 
in the wounded nucleon model. The difference between 
the exact Monte Carlo results and the analytic formulas 
is due to the presence of correlations between the location 
of sources in Glauber- like models. These correlations re- 
sult from a rather simple mechanism mentioned at the 
beginning of Sect. I VIII a nucleon from nucleus A may 
wound several nucleons from nucleus B. This results in 
some clustering, hence correlations, of the locations of 
the wounded nucleons. In the derivation of the analytic 
formulas we have resorted to the central limit theorem, 
hence all correlations were neglected. The difference be- 
tween the full and uncorrelated analytic results in Fig. [15] 
display the significance of the correlations. 

The correlations between the locations of sources re- 
sult in a decrease of the effective number of sources n, 
so the full result for e* (the monotonically rising curves) 
is naturally above the uncorrelated analytic result. The 
behavior of the two curves is similar and the relative dif- 
ference is at the level of 10-15%. For the scaled standard 
deviation, Ae*/e*, the comparison is more complicated. 
At low values of b the two curves are very close, at in- 
termediate b the calculation with correlations is higher, 
while at peripheric b it is lower than the uncorrelated 
case. We conclude that at central collisions Ae*/e* is 
not sensitive to correlations. 




10 12 



b [fm] 



FIG. 15: (Color online) Comparison of the Monte Carlo cal- 
culation of e* = £2,2 (the rising curves) and Ae*/e* in the 
wounded nucleon model for gold-gold collisions (solid lines) 
and the analytic formulas (|D10|) with / = 2 (dashed lines). 
The analytic formulas neglect correlations between the loca- 
tion of sources present in the full calculation. 



In Monte Carlo simulations the correlations may be 
artificially removed by taking a very large cross section 
cr w , in which case all the nucleons get wounded and the 
correlations between the locations of sources disappear. 
This may be used for testing purposes. In that case the 
two calculations of Fig. [TBI overlap. 

An analytic inclusion of correlations into the frame- 
work based on the central limit theorem is difficult and 
it is more productive to simply perform the simulations. 
However, the analytic formulas (|D10[) bare significance 
not only at the formal level, which helps to understand 
the nature of the chosen statistical measures. There may 
be some models where the correlations are largely re- 
duced compared to the wounded nucleon model, or ab- 
sent. Then the evaluation of e* and its variance are sim- 
ply made by computing the moments ^2,0 , h,o, h,2, and 
^4,4 of the fixed- axes distribution and carrying out a trun- 
cated series in Eq. ( plO ). We note that, amusingly, the 
CGC calculation of As* /e* shown in Fig. ([7]) agrees sur- 
prisingly well with the uncorrelated result from Fig. (|15p . 
This hints that the CGC approach of Ref. [l9[ has un- 
correlated sources. 

Next, we derive expressions for the variable-axes pro- 
files fi(p) in the absence of particle correlations. These 
profiles correspond to inclusive distributions uninte- 
grated over the p variable of a selected particle, namely 



27r PJ7(p) = J dej) j dxi . . . dx n f{xi) . . . f(x n ) x 

n 

]T S( Pm - p)S(cb m - <f>) cos[l(<f> - n (Dll) 



m— 1 



where the single-particle distributions f(xi) are normal- 
ized to unity. The inclusive distribution is normalized to 
n, hence Iq.q = J 27rp/g {p)dp = n. Since all particles 
have equal distributions f(xi), we may relabel particles 
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setting, for instance, x n = (p, <fi) and rewrite Eq. (|D1 1[) as 
I 



flip) 



7rf(p><t>) [ dxi...dx n -if(xi)...f(x n -i)- 
2ir J q 



cos(l(j>) ^2 Pm cos(Z0 m ) + shx(l<f>) ^2 Pm sin(^ri 



(D12) 



where we have used the definitions (|7.3p and (|D5|) . Similarly to the analysis of the moments of Sect. IVIIi for sufficiently 
large values of n the variables 



1 n-l 1 l n-l 1 

Yi = - p k cos(ld>j) =Y t p k cos</> = q' cos(a), X' = - p h - sin(/0,) = X; p k sin</> = </sin(a), 

77. — ' -' 77. 77 — ' 77. 



(D13) 



follow the normal distribution with rescaled parameters Y( = 1L ^-Yi and ct^? x = 1L ^-&y l X[ . Therefore 



dq'da q' 
27ra' Y a' X[ q 

dq'da 



exp 



Y/ 2 g'Y~zcosa 



2a' 



5'q' 2 sin 2 a 



[//' + nq' cos(/</> - a)] 



2-no' y o' Xi y/q'2 + 2q' p k cos(l</> -a)/n + p 2k /r 



exp 



q' 2 + Yi q'Yi cos a 



2a' 



5'q' sin 2 a 



[p fc + nq' cos(l(j) — a)] 



For the case of central collisions, where 5' = 0, Yi = 0, and a'y l = a' x l , formula (|D14p simplifies into 



flip) = hip) 

= Up) 
= Hp) 



dq'd/3 



2-Ka'y v <?' 2 + 2q'p k /n cos /3 + p 2k /n- 



: exp 



27rcr^ 



27rcr: 



/ 2 



■ exp 



exp 



/2 




q 


2n 







(1 + nq/p k ) E 



Anqp k 
(p k + nq) 



[p k + nq' cos 0\ 
r ) +(l-nq/p k )K 



(D14) 



(D15) 



Anqp k 
(p k + nqf 



2a 



/2 



■2k- 



8n 2 q 



72 +■••). 



where £ and if denote the elliptic integrals of the second profile /q. 



and third kind. Since q' is of the order of a'y l ~ 1/y/n, 
subsequent terms in the expansion denoted by . . . are 



Also note, that the power in the multiplying factor 
p k simply reflects the (arbitrarily) chosen power for the 



suppressed with powers of n. Hence, in the large-n limit averaging; thus is a matter of methodology (c/. Tabic HJ) 
we may retain only the first term in the expansion (the 



unity). The same result is obtained by first expanding 
q'/q in inverse powers of n and then carrying the inte- 
gration over (3. The remaining integral over q' is trivial, 
yielding the final expression for the central case in the 
absence of correlations: 



We now return to the non-central case of Eq. (|D14|> . 
We first expand the following piece in the inverse powers, 
of n, retaining the first two: 

(p k + nq' cos(l(/) — a))q' 



2 V nhkA 



■p k fo(p), (6 = 0). (D16) 



\/q' 2 + 2q' p k /ncos(lcj) — a) + p 2k /n 2 

= nq' cos(Z0 - a) + p k sin 2 (Z0 - a) + . . . (D17) 



Remarkably, in this case all variable-axes profiles are 

equal to one another and depend only on the monopole We may then carry the integration over <j>, which gives 
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flip) 



dq'da 



2ira' 



■ exp 



?' 2 + 5f 



q'Yi cos a 



2a' 2 



c/ /2 • 2 

o q sin a 



Yi 



Y, 



nq'fi(p) cos(a) + -p k f (p) - ^/M/>) cos(2a) 



(D18) 



We note the presence of three fixed-axes profiles: fi(p), 
p k fo(p), and p k f2i(p)- The integration over a may be 
done similarly to the case of the moments, via expan- 
sion in the Bessel functions and then carrying out the q' 



integration. The result, involving various confluent hy- 
pergeometric functions, is rather lengthy hence we do not 
list it here. 
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